import pandas as pd
import numpy as np
import random
import matplotlib.pyplot as plt
import os
import time

# 定义模拟退火算法类
class SimulatedAnnealing:
    def __init__(self, initial_params, param_ranges, max_iterations, initial_temp, cooling_rate, data):
        self.params = initial_params  # 初始参数 [k, a, b]
        self.param_ranges = param_ranges  # 参数范围
        self.max_iterations = max_iterations  # 最大迭代次数
        self.temp = initial_temp  # 初始温度
        self.cooling_rate = cooling_rate  # 降温速率
        self.data = data  # 数据
        self.best_solution = initial_params  # 最优解
        self.best_error, self.best_mae, self.best_rmse = self.objective_function(self.params)  # 最优误差
        self.error_history = []  # 记录每次迭代的误差

    # 目标函数
    def objective_function(self, params):
        k, a, b = params
        P_actual = self.data['P']
        f = self.data['f']
        B = self.data['B']
        P_predicted = k * (f ** a) * (B ** b)

        # 使用相对误差
        relative_error = np.mean(((P_actual - P_predicted) / P_actual) ** 2)
        # 平均绝对误差 (MAE)
        mae = np.mean(np.abs(P_actual - P_predicted))
        # 均方根误差 (RMSE)
        rmse = np.sqrt(np.mean((P_actual - P_predicted) ** 2))

        return relative_error, mae, rmse

    # 执行一次模拟退火迭代
    def run(self):
        for iteration in range(self.max_iterations):
            # 生成新的解
            new_params = self.generate_new_solution(self.params)
            new_error, new_mae, new_rmse = self.objective_function(new_params)

            # 比较新解与当前解
            if new_error < self.best_error:
                self.params = new_params
                self.best_error = new_error
                self.best_mae = new_mae
                self.best_rmse = new_rmse
                self.best_solution = new_params
            else:
                # 按一定概率接受较差解
                if self.accept_solution(new_error):
                    self.params = new_params

            # 降低温度
            self.temp *= self.cooling_rate
            self.error_history.append(self.best_error)

        return self.best_solution, self.best_error, self.best_mae, self.best_rmse

    # 生成新解
    def generate_new_solution(self, current_params):
        new_params = []
        for i, (low, high) in enumerate(self.param_ranges):
            perturbation = random.uniform(-0.1, 0.1)  # 在参数附近进行扰动
            new_param = current_params[i] + perturbation
            # 确保新参数在指定范围内
            new_param = max(low, min(high, new_param))
            new_params.append(new_param)
        return new_params

    # 接受新解的概率
    def accept_solution(self, new_error):
        delta_error = new_error - self.best_error
        acceptance_probability = np.exp(-delta_error / self.temp)
        return random.random() < acceptance_probability

    # 绘制迭代误差曲线
    def plot_error_history(self, output_dir, run):
        plt.figure()
        plt.plot(self.error_history, label='errorline')
        plt.title(f'num_error (Num: {run})')
        plt.xlabel('num')
        plt.ylabel('error')
        plt.legend()
        plt.grid()
        plt.savefig(os.path.join(output_dir, f'num_error{run}.png'))
        plt.close()

    # 绘制拟合曲线
    def plot_fit_curve(self, output_dir, run):
        k, a, b = self.best_solution
        f = self.data['f']
        B = self.data['B']
        P_actual = self.data['P']
        P_predicted = k * (f ** a) * (B ** b)

        plt.figure(figsize=(10, 6))
        plt.scatter(P_predicted, P_actual, label='P', alpha=0.5)
        plt.plot([min(P_actual), max(P_actual)], [min(P_actual), max(P_actual)], 'r--', label='Y=X', linewidth=2)
        plt.title(f'sim_num (run_num: {run})')
        plt.xlabel('preP')
        plt.ylabel('trueP')
        plt.legend()
        plt.grid()
        plt.savefig(os.path.join(output_dir, f'trueP VS preP{run}.png'))
        plt.close()


# 读取CSV文件路径
file_paths = [
    r'D:\2024华为杯数学建模\2024C题-数据驱动与电磁损耗\C题\附件一数据拆解\Q2\Q2磁芯材料1正弦波形.csv',
]

# 合并所有数据
data_frames = []
for file_path in file_paths:
    if not os.path.exists(file_path):
        raise FileNotFoundError(f"文件不存在: {file_path}")
    try:
        df = pd.read_csv(file_path, skiprows=1, encoding='gbk')
    except UnicodeDecodeError:
        df = pd.read_csv(file_path, skiprows=1, encoding='ISO-8859-1')
    data_frames.append(df)

data = pd.concat(data_frames, ignore_index=True)
num_columns = data.shape[1]
new_column_names = [f'Column_{i+1}' for i in range(num_columns)]
data.columns = new_column_names
data.rename(columns={'Column_2': 'f', 'Column_3': 'P', 'Column_5': 'B'}, inplace=True)

# 定义参数范围
param_ranges = [
    (0.1, 10),  # k 的范围
    (1, 3),     # a 的范围
    (2, 3)      # b 的范围
]

# 记录运行时间
start_time = time.time()

# 输出目录
output_dir = r'D:\2024华为杯数学建模\2024C题-数据驱动与电磁损耗\C_code\Q2_code\Q2可视化结果\Q2参数拟合图示'
os.makedirs(output_dir, exist_ok=True)  # 创建目录

# 运行模拟退火算法 50次
results = []
for run in range(50):
    # 初始化参数
    initial_params = [
        random.uniform(*param_ranges[0]),  # k 的初始值
        random.uniform(*param_ranges[1]),  # a 的初始值
        random.uniform(*param_ranges[2])   # b 的初始值
    ]

    # 实例化模拟退火算法
    sa = SimulatedAnnealing(
        initial_params=initial_params,
        param_ranges=param_ranges,
        max_iterations=1000,  # 设置迭代次数
        initial_temp=1000,    # 初始温度
        cooling_rate=0.95,    # 降温速率
        data=data
    )

    # 运行模拟退火算法
    best_params, best_error, best_mae, best_rmse = sa.run()

    # 保存结果
    results.append((run + 1, best_params[0], best_params[1], best_params[2], best_error, best_mae, best_rmse))

    # 绘制迭代误差曲线和拟合曲线
    sa.plot_error_history(output_dir, run + 1)
    sa.plot_fit_curve(output_dir, run + 1)

end_time = time.time()

# 输出运行时间
print(f"运行时间: {end_time - start_time:.2f} 秒")

# 保存最优结果到CSV文件
output_csv_path = os.path.join(output_dir, '最优拟合结果多次.csv')
result_df = pd.DataFrame(results, columns=['运行次数', 'k', 'a', 'b', '最小误差', 'MAE', 'RMSE'])
result_df.to_csv(output_csv_path, index=False)

print(f"所有结果已保存到: {output_csv_path}")